home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / STUTTGART / LANG / C / GCC / V2-4-5 / GPPLIBSR00 / cc / floatconv < prev    next >
Text File  |  1993-12-08  |  76KB  |  2,405 lines
  1. #include "ioprivate.h"
  2. #ifdef USE_DTOA
  3. /****************************************************************
  4.  *
  5.  * The author of this software is David M. Gay.
  6.  *
  7.  * Copyright (c) 1991 by AT&T.
  8.  *
  9.  * Permission to use, copy, modify, and distribute this software for any
  10.  * purpose without fee is hereby granted, provided that this entire notice
  11.  * is included in all copies of any software which is or includes a copy
  12.  * or modification of this software and in all copies of the supporting
  13.  * documentation for such software.
  14.  *
  15.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  16.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  17.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  18.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  19.  *
  20.  ***************************************************************/
  21.  
  22. /* Please send bug reports to
  23.         David M. Gay
  24.         AT&T Bell Laboratories, Room 2C-463
  25.         600 Mountain Avenue
  26.         Murray Hill, NJ 07974-2070
  27.         U.S.A.
  28.         dmg@research.att.com or research!dmg
  29.  */
  30.  
  31. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  32.  *
  33.  * This strtod returns a nearest machine number to the input decimal
  34.  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  35.  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  36.  * biased rounding (add half and chop).
  37.  *
  38.  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  39.  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  40.  *
  41.  * Modifications:
  42.  *
  43.  *      1. We only require IEEE, IBM, or VAX double-precision
  44.  *              arithmetic (not IEEE double-extended).
  45.  *      2. We get by with floating-point arithmetic in a case that
  46.  *              Clinger missed -- when we're computing d * 10^n
  47.  *              for a small integer d and the integer n is not too
  48.  *              much larger than 22 (the maximum integer k for which
  49.  *              we can represent 10^k exactly), we may be able to
  50.  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  51.  *      3. Rather than a bit-at-a-time adjustment of the binary
  52.  *              result in the hard case, we use floating-point
  53.  *              arithmetic to determine the adjustment to within
  54.  *              one bit; only in really hard cases do we need to
  55.  *              compute a second residual.
  56.  *      4. Because of 3., we don't need a large table of powers of 10
  57.  *              for ten-to-e (just some small tables, e.g. of 10^k
  58.  *              for 0 <= k <= 22).
  59.  */
  60.  
  61. /*
  62.  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  63.  *      significant byte has the lowest address.
  64.  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  65.  *      significant byte has the lowest address.
  66.  * #define Sudden_Underflow for IEEE-format machines without gradual
  67.  *      underflow (i.e., that flush to zero on underflow).
  68.  * #define IBM for IBM mainframe-style floating-point arithmetic.
  69.  * #define VAX for VAX-style floating-point arithmetic.
  70.  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
  71.  * #define No_leftright to omit left-right logic in fast floating-point
  72.  *      computation of dtoa.
  73.  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  74.  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  75.  *      that use extended-precision instructions to compute rounded
  76.  *      products and quotients) with IBM.
  77.  * #define ROUND_BIASED for IEEE-format with biased rounding.
  78.  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  79.  *      products but inaccurate quotients, e.g., for Intel i860.
  80.  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
  81.  *      integer arithmetic.  Whether this speeds things up or slows things
  82.  *      down depends on the machine and the number being converted.
  83.  * #define KR_headers for old-style C function headers.
  84.  */
  85.  
  86. #ifdef DEBUG
  87. extern "C"
  88.   {
  89.   #include <stdio.h>
  90.   }
  91.  
  92. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  93. #endif
  94.  
  95. extern "C"
  96.   {
  97.   #include <stdlib.h>
  98.   #include <string.h>
  99.   }
  100.  
  101. #define CONST const
  102.  
  103. extern "C"
  104.   {
  105.   #include <errno.h>
  106.   #include <float.h>
  107.   #ifndef __MATH_H__
  108.   #include <math.h>
  109.   #endif
  110.   }
  111.  
  112.  
  113. #ifdef Unsigned_Shifts
  114. #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
  115. #else
  116. #define Sign_Extend(a,b) /*no-op*/
  117. #endif
  118.  
  119. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  120.  
  121. #if FLT_RADIX==16
  122. #define IBM
  123. #elif DBL_MANT_DIG==56
  124. #define VAX
  125. #elif DBL_MANT_DIG==53 && DBL_MAX_10_EXP==308
  126. #define IEEE_Unknown
  127. #else
  128. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  129. #endif
  130. #endif
  131.  
  132. #ifdef IEEE_8087
  133. #define HIWORD 1
  134. #define LOWORD 0
  135. #define TEST_ENDIANNESS  /* nothing */
  136. #elif defined(IEEE_MC68k)
  137. #define HIWORD 0
  138. #define LOWORD 1
  139. #define TEST_ENDIANNESS  /* nothing */
  140. #else
  141. static int HIWORD = -1, LOWORD;
  142. static void test_endianness()
  143. {
  144.     union doubleword {
  145.     double d;
  146.     unsigned long u[2];
  147.     } dw;
  148.     dw.d = 10;
  149.     if (dw.u[0] != 0) /* big-endian */
  150.     HIWORD=0, LOWORD=1;
  151.     else
  152.     HIWORD=1, LOWORD=0;
  153. }
  154. #define TEST_ENDIANNESS  if (HIWORD<0) test_endianness();
  155. #endif
  156. #define word0(x) ((unsigned long *)&x)[HIWORD]
  157. #define word1(x) ((unsigned long *)&x)[LOWORD]
  158.  
  159. /* The following definition of Storeinc is appropriate for MIPS processors. */
  160. #if defined(IEEE_8087) + defined(VAX)
  161. #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
  162. ((unsigned short *)a)[0] = (unsigned short)c, a++)
  163. #elif defined(IEEE_MC68k)
  164. #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
  165. ((unsigned short *)a)[1] = (unsigned short)c, a++)
  166. #else
  167. #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  168. #endif
  169.  
  170. /* #define P DBL_MANT_DIG */
  171. /* Ten_pmax = floor(P*log(2)/log(5)) */
  172. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  173. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  174. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  175.  
  176. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_Unknown)
  177. #define Exp_shift  20
  178. #define Exp_shift1 20
  179. #define Exp_msk1    0x100000
  180. #define Exp_msk11   0x100000
  181. #define Exp_mask  0x7ff00000
  182. #define P 53
  183. #define Bias 1023
  184. #define IEEE_Arith
  185. #define Emin (-1022)
  186. #define Exp_1  0x3ff00000
  187. #define Exp_11 0x3ff00000
  188. #define Ebits 11
  189. #define Frac_mask  0xfffff
  190. #define Frac_mask1 0xfffff
  191. #define Ten_pmax 22
  192. #define Bletch 0x10
  193. #define Bndry_mask  0xfffff
  194. #define Bndry_mask1 0xfffff
  195. #define LSB 1
  196. #define Sign_bit 0x80000000
  197. #define Log2P 1
  198. #define Tiny0 0
  199. #define Tiny1 1
  200. #define Quick_max 14
  201. #define Int_max 14
  202. #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
  203. #else
  204. #undef  Sudden_Underflow
  205. #define Sudden_Underflow
  206. #ifdef IBM
  207. #define Exp_shift  24
  208. #define Exp_shift1 24
  209. #define Exp_msk1   0x1000000
  210. #define Exp_msk11  0x1000000
  211. #define Exp_mask  0x7f000000
  212. #define P 14
  213. #define Bias 65
  214. #define Exp_1  0x41000000
  215. #define Exp_11 0x41000000
  216. #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
  217. #define Frac_mask  0xffffff
  218. #define Frac_mask1 0xffffff
  219. #define Bletch 4
  220. #define Ten_pmax 22
  221. #define Bndry_mask  0xefffff
  222. #define Bndry_mask1 0xffffff
  223. #define LSB 1
  224. #define Sign_bit 0x80000000
  225. #define Log2P 4
  226. #define Tiny0 0x100000
  227. #define Tiny1 0
  228. #define Quick_max 14
  229. #define Int_max 15
  230. #else /* VAX */
  231. #define Exp_shift  23
  232. #define Exp_shift1 7
  233. #define Exp_msk1    0x80
  234. #define Exp_msk11   0x800000
  235. #define Exp_mask  0x7f80
  236. #define P 56
  237. #define Bias 129
  238. #define Exp_1  0x40800000
  239. #define Exp_11 0x4080
  240. #define Ebits 8
  241. #define Frac_mask  0x7fffff
  242. #define Frac_mask1 0xffff007f
  243. #define Ten_pmax 24
  244. #define Bletch 2
  245. #define Bndry_mask  0xffff007f
  246. #define Bndry_mask1 0xffff007f
  247. #define LSB 0x10000
  248. #define Sign_bit 0x8000
  249. #define Log2P 1
  250. #define Tiny0 0x80
  251. #define Tiny1 0
  252. #define Quick_max 15
  253. #define Int_max 15
  254. #endif
  255. #endif
  256.  
  257. #ifndef IEEE_Arith
  258. #define ROUND_BIASED
  259. #endif
  260.  
  261. #ifdef RND_PRODQUOT
  262. #define rounded_product(a,b) a = rnd_prod(a, b)
  263. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  264. extern double rnd_prod(double, double), rnd_quot(double, double);
  265. #else
  266. #define rounded_product(a,b) a *= b
  267. #define rounded_quotient(a,b) a /= b
  268. #endif
  269.  
  270. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  271. #define Big1 0xffffffff
  272.  
  273. #ifndef Just_16
  274. /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
  275.  * This makes some inner loops simpler and sometimes saves work
  276.  * during multiplications, but it often seems to make things slightly
  277.  * slower.  Hence the default is now to store 32 bits per long.
  278.  */
  279. #ifndef Pack_32
  280. #define Pack_32
  281. #endif
  282. #endif
  283.  
  284. #define Kmax 15
  285.  
  286. extern "C" double _Xstrtod(const char *s00, char **se);
  287. extern "C" char *dtoa(double d, int mode, int ndigits,
  288.                         int *decpt, int *sign, char **rve);
  289.  
  290.  struct
  291. Bigint {
  292.         struct Bigint *next;
  293.         int k, maxwds, sign, wds;
  294.         unsigned long x[1];
  295.         };
  296.  
  297.  typedef struct Bigint Bigint;
  298.  
  299.  static Bigint *freelist[Kmax+1];
  300.  
  301.  static Bigint *
  302. Balloc
  303. #ifdef KR_headers
  304.         (k) int k;
  305. #else
  306.         (int k)
  307. #endif
  308. {
  309.         int x;
  310.         Bigint *rv;
  311.  
  312.         if (rv = freelist[k]) {
  313.                 freelist[k] = rv->next;
  314.                 }
  315.         else {
  316.                 x = 1 << k;
  317.                 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
  318.                 rv->k = k;
  319.                 rv->maxwds = x;
  320.                 }
  321.         rv->sign = rv->wds = 0;
  322.         return rv;
  323.         }
  324.  
  325.  static void
  326. Bfree
  327. #ifdef KR_headers
  328.         (v) Bigint *v;
  329. #else
  330.         (Bigint *v)
  331. #endif
  332. {
  333.         if (v) {
  334.                 v->next = freelist[v->k];
  335.                 freelist[v->k] = v;
  336.                 }
  337.         }
  338.  
  339. #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
  340. y->wds*sizeof(long) + 2*sizeof(int))
  341.  
  342.  static Bigint *
  343. multadd
  344. #ifdef KR_headers
  345.         (b, m, a) Bigint *b; int m, a;
  346. #else
  347.         (Bigint *b, int m, int a)       /* multiply by m and add a */
  348. #endif
  349. {
  350.         int i, wds;
  351.         unsigned long *x, y;
  352. #ifdef Pack_32
  353.         unsigned long xi, z;
  354. #endif
  355.         Bigint *b1;
  356.  
  357.         wds = b->wds;
  358.         x = b->x;
  359.         i = 0;
  360.         do {
  361. #ifdef Pack_32
  362.                 xi = *x;
  363.                 y = (xi & 0xffff) * m + a;
  364.                 z = (xi >> 16) * m + (y >> 16);
  365.                 a = (int)(z >> 16);
  366.                 *x++ = (z << 16) + (y & 0xffff);
  367. #else
  368.                 y = *x * m + a;
  369.                 a = (int)(y >> 16);
  370.                 *x++ = y & 0xffff;
  371. #endif
  372.                 }
  373.                 while(++i < wds);
  374.         if (a) {
  375.                 if (wds >= b->maxwds) {
  376.                         b1 = Balloc(b->k+1);
  377.                         Bcopy(b1, b);
  378.                         Bfree(b);
  379.                         b = b1;
  380.                         }
  381.                 b->x[wds++] = a;
  382.                 b->wds = wds;
  383.                 }
  384.         return b;
  385.         }
  386.  
  387.  static Bigint *
  388. s2b
  389. #ifdef KR_headers
  390.         (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
  391. #else
  392.         (CONST char *s, int nd0, int nd, unsigned long y9)
  393. #endif
  394. {
  395.         Bigint *b;
  396.         int i, k;
  397.         long x, y;
  398.  
  399.         x = (nd + 8) / 9;
  400.         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
  401. #ifdef Pack_32
  402.         b = Balloc(k);
  403.         b->x[0] = y9;
  404.         b->wds = 1;
  405. #else
  406.         b = Balloc(k+1);
  407.         b->x[0] = y9 & 0xffff;
  408.         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
  409. #endif
  410.  
  411.         i = 9;
  412.         if (9 < nd0) {
  413.                 s += 9;
  414.                 do b = multadd(b, 10, *s++ - '0');
  415.                         while(++i < nd0);
  416.                 s++;
  417.                 }
  418.         else
  419.                 s += 10;
  420.         for(; i < nd; i++)
  421.                 b = multadd(b, 10, *s++ - '0');
  422.         return b;
  423.         }
  424.  
  425.  static int
  426. hi0bits
  427. #ifdef KR_headers
  428.         (x) register unsigned long x;
  429. #else
  430.         (register unsigned long x)
  431. #endif
  432. {
  433.         register int k = 0;
  434.  
  435.         if (!(x & 0xffff0000)) {
  436.                 k = 16;
  437.                 x <<= 16;
  438.                 }
  439.         if (!(x & 0xff000000)) {
  440.                 k += 8;
  441.                 x <<= 8;
  442.                 }
  443.         if (!(x & 0xf0000000)) {
  444.                 k += 4;
  445.                 x <<= 4;
  446.                 }
  447.         if (!(x & 0xc0000000)) {
  448.                 k += 2;
  449.                 x <<= 2;
  450.                 }
  451.         if (!(x & 0x80000000)) {
  452.                 k++;
  453.                 if (!(x & 0x40000000))
  454.                         return 32;
  455.                 }
  456.         return k;
  457.         }
  458.  
  459.  static int
  460. lo0bits
  461. #ifdef KR_headers
  462.         (y) unsigned long *y;
  463. #else
  464.         (unsigned long *y)
  465. #endif
  466. {
  467.         register int k;
  468.         register unsigned long x = *y;
  469.  
  470.         if (x & 7) {
  471.                 if (x & 1)
  472.                         return 0;
  473.                 if (x & 2) {
  474.                         *y = x >> 1;
  475.                         return 1;
  476.                         }
  477.                 *y = x >> 2;
  478.                 return 2;
  479.                 }
  480.         k = 0;
  481.         if (!(x & 0xffff)) {
  482.                 k = 16;
  483.                 x >>= 16;
  484.                 }
  485.         if (!(x & 0xff)) {
  486.                 k += 8;
  487.                 x >>= 8;
  488.                 }
  489.         if (!(x & 0xf)) {
  490.                 k += 4;
  491.                 x >>= 4;
  492.                 }
  493.         if (!(x & 0x3)) {
  494.                 k += 2;
  495.                 x >>= 2;
  496.                 }
  497.         if (!(x & 1)) {
  498.                 k++;
  499.                 x >>= 1;
  500.                 if (!x & 1)
  501.                         return 32;
  502.                 }
  503.         *y = x;
  504.         return k;
  505.         }
  506.  
  507.  static Bigint *
  508. i2b
  509. #ifdef KR_headers
  510.         (i) int i;
  511. #else
  512.         (int i)
  513. #endif
  514. {
  515.         Bigint *b;
  516.  
  517.         b = Balloc(1);
  518.         b->x[0] = i;
  519.         b->wds = 1;
  520.         return b;
  521.         }
  522.  
  523.  static Bigint *
  524. mult
  525. #ifdef KR_headers
  526.         (a, b) Bigint *a, *b;
  527. #else
  528.         (Bigint *a, Bigint *b)
  529. #endif
  530. {
  531.         Bigint *c;
  532.         int k, wa, wb, wc;
  533.         unsigned long carry, y, z;
  534.         unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  535. #ifdef Pack_32
  536.         unsigned long z2;
  537. #endif
  538.  
  539.         if (a->wds < b->wds) {
  540.                 c = a;
  541.                 a = b;
  542.                 b = c;
  543.                 }
  544.         k = a->k;
  545.         wa = a->wds;
  546.         wb = b->wds;
  547.         wc = wa + wb;
  548.         if (wc > a->maxwds)
  549.                 k++;
  550.         c = Balloc(k);
  551.         for(x = c->x, xa = x + wc; x < xa; x++)
  552.                 *x = 0;
  553.         xa = a->x;
  554.         xae = xa + wa;
  555.         xb = b->x;
  556.         xbe = xb + wb;
  557.         xc0 = c->x;
  558. #ifdef Pack_32
  559.         for(; xb < xbe; xb++, xc0++) {
  560.                 if (y = *xb & 0xffff) {
  561.                         x = xa;
  562.                         xc = xc0;
  563.                         carry = 0;
  564.                         do {
  565.                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  566.                                 carry = z >> 16;
  567.                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  568.                                 carry = z2 >> 16;
  569.                                 Storeinc(xc, z2, z);
  570.                                 }
  571.                                 while(x < xae);
  572.                         *xc = carry;
  573.                         }
  574.                 if (y = *xb >> 16) {
  575.                         x = xa;
  576.                         xc = xc0;
  577.                         carry = 0;
  578.                         z2 = *xc;
  579.                         do {
  580.                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  581.                                 carry = z >> 16;
  582.                                 Storeinc(xc, z, z2);
  583.                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  584.                                 carry = z2 >> 16;
  585.                                 }
  586.                                 while(x < xae);
  587.                         *xc = z2;
  588.                         }
  589.                 }
  590. #else
  591.         for(; xb < xbe; xc0++) {
  592.                 if (y = *xb++) {
  593.                         x = xa;
  594.                         xc = xc0;
  595.                         carry = 0;
  596.                         do {
  597.                                 z = *x++ * y + *xc + carry;
  598.                                 carry = z >> 16;
  599.                                 *xc++ = z & 0xffff;
  600.                                 }
  601.                                 while(x < xae);
  602.                         *xc = carry;
  603.                         }
  604.                 }
  605. #endif
  606.         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  607.         c->wds = wc;
  608.         return c;
  609.         }
  610.  
  611.  static Bigint *p5s;
  612.  
  613.  static Bigint *
  614. pow5mult
  615. #ifdef KR_headers
  616.         (b, k) Bigint *b; int k;
  617. #else
  618.         (Bigint *b, int k)
  619. #endif
  620. {
  621.         Bigint *b1, *p5, *p51;
  622.         int i;
  623.         static int p05[3] = { 5, 25, 125 };
  624.  
  625.         if (i = k & 3)
  626.                 b = multadd(b, p05[i-1], 0);
  627.  
  628.         if (!(k >>= 2))
  629.                 return b;
  630.         if (!(p5 = p5s)) {
  631.                 /* first time */
  632.                 p5 = p5s = i2b(625);
  633.                 p5->next = 0;
  634.                 }
  635.         for(;;) {
  636.                 if (k & 1) {
  637.                         b1 = mult(b, p5);
  638.                         Bfree(b);
  639.                         b = b1;
  640.                         }
  641.                 if (!(k >>= 1))
  642.                         break;
  643.                 if (!(p51 = p5->next)) {
  644.                         p51 = p5->next = mult(p5,p5);
  645.                         p51->next = 0;
  646.                         }
  647.                 p5 = p51;
  648.                 }
  649.         return b;
  650.         }
  651.  
  652.  static Bigint *
  653. lshift
  654. #ifdef KR_headers
  655.         (b, k) Bigint *b; int k;
  656. #else
  657.         (Bigint *b, int k)
  658. #endif
  659. {
  660.         int i, k1, n, n1;
  661.         Bigint *b1;
  662.         unsigned long *x, *x1, *xe, z;
  663.  
  664. #ifdef Pack_32
  665.         n = k >> 5;
  666. #else
  667.         n = k >> 4;
  668. #endif
  669.         k1 = b->k;
  670.         n1 = n + b->wds + 1;
  671.         for(i = b->maxwds; n1 > i; i <<= 1)
  672.                 k1++;
  673.         b1 = Balloc(k1);
  674.         x1 = b1->x;
  675.         for(i = 0; i < n; i++)
  676.                 *x1++ = 0;
  677.         x = b->x;
  678.         xe = x + b->wds;
  679. #ifdef Pack_32
  680.         if (k &= 0x1f) {
  681.                 k1 = 32 - k;
  682.                 z = 0;
  683.                 do {
  684.                         *x1++ = *x << k | z;
  685.                         z = *x++ >> k1;
  686.                         }
  687.                         while(x < xe);
  688.                 if (*x1 = z)
  689.                         ++n1;
  690.                 }
  691. #else
  692.         if (k &= 0xf) {
  693.                 k1 = 16 - k;
  694.                 z = 0;
  695.                 do {
  696.                         *x1++ = *x << k  & 0xffff | z;
  697.                         z = *x++ >> k1;
  698.                         }
  699.                         while(x < xe);
  700.                 if (*x1 = z)
  701.                         ++n1;
  702.                 }
  703. #endif
  704.         else do
  705.                 *x1++ = *x++;
  706.                 while(x < xe);
  707.         b1->wds = n1 - 1;
  708.         Bfree(b);
  709.         return b1;
  710.         }
  711.  
  712.  static int
  713. cmp
  714. #ifdef KR_headers
  715.         (a, b) Bigint *a, *b;
  716. #else
  717.         (Bigint *a, Bigint *b)
  718. #endif
  719. {
  720.         unsigned long *xa, *xa0, *xb, *xb0;
  721.         int i, j;
  722.  
  723.         i = a->wds;
  724.         j = b->wds;
  725. #ifdef DEBUG
  726.         if (i > 1 && !a->x[i-1])
  727.                 Bug("cmp called with a->x[a->wds-1] == 0");
  728.         if (j > 1 && !b->x[j-1])
  729.                 Bug("cmp called with b->x[b->wds-1] == 0");
  730. #endif
  731.         if (i -= j)
  732.                 return i;
  733.         xa0 = a->x;
  734.         xa = xa0 + j;
  735.         xb0 = b->x;
  736.         xb = xb0 + j;
  737.         for(;;) {
  738.                 if (*--xa != *--xb)
  739.                         return *xa < *xb ? -1 : 1;
  740.                 if (xa <= xa0)
  741.                         break;
  742.                 }
  743.         return 0;
  744.         }
  745.  
  746.  static Bigint *
  747. diff
  748. #ifdef KR_headers
  749.         (a, b) Bigint *a, *b;
  750. #else
  751.         (Bigint *a, Bigint *b)
  752. #endif
  753. {
  754.         Bigint *c;
  755.         int i, wa, wb;
  756.         long borrow, y; /* We need signed shifts here. */
  757.         unsigned long *xa, *xae, *xb, *xbe, *xc;
  758. #ifdef Pack_32
  759.         long z;
  760. #endif
  761.  
  762.         i = cmp(a,b);
  763.         if (!i) {
  764.                 c = Balloc(0);
  765.                 c->wds = 1;
  766.                 c->x[0] = 0;
  767.                 return c;
  768.                 }
  769.         if (i < 0) {
  770.                 c = a;
  771.                 a = b;
  772.                 b = c;
  773.                 i = 1;
  774.                 }
  775.         else
  776.                 i = 0;
  777.         c = Balloc(a->k);
  778.         c->sign = i;
  779.         wa = a->wds;
  780.         xa = a->x;
  781.         xae = xa + wa;
  782.         wb = b->wds;
  783.         xb = b->x;
  784.         xbe = xb + wb;
  785.         xc = c->x;
  786.         borrow = 0;
  787. #ifdef Pack_32
  788.         do {
  789.                 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
  790.                 borrow = y >> 16;
  791.                 Sign_Extend(borrow, y);
  792.                 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
  793.                 borrow = z >> 16;
  794.                 Sign_Extend(borrow, z);
  795.                 Storeinc(xc, z, y);
  796.                 }
  797.                 while(xb < xbe);
  798.         while(xa < xae) {
  799.                 y = (*xa & 0xffff) + borrow;
  800.                 borrow = y >> 16;
  801.                 Sign_Extend(borrow, y);
  802.                 z = (*xa++ >> 16) + borrow;
  803.                 borrow = z >> 16;
  804.                 Sign_Extend(borrow, z);
  805.                 Storeinc(xc, z, y);
  806.                 }
  807. #else
  808.         do {
  809.                 y = *xa++ - *xb++ + borrow;
  810.                 borrow = y >> 16;
  811.                 Sign_Extend(borrow, y);
  812.                 *xc++ = y & 0xffff;
  813.                 }
  814.                 while(xb < xbe);
  815.         while(xa < xae) {
  816.                 y = *xa++ + borrow;
  817.                 borrow = y >> 16;
  818.                 Sign_Extend(borrow, y);
  819.                 *xc++ = y & 0xffff;
  820.                 }
  821. #endif
  822.         while(!*--xc)
  823.                 wa--;
  824.         c->wds = wa;
  825.         return c;
  826.         }
  827.  
  828.  static double
  829. ulp
  830. #ifdef KR_headers
  831.         (x) double x;
  832. #else
  833.         (double x)
  834. #endif
  835. {
  836.         register long L;
  837.         double a;
  838.  
  839.         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
  840. #ifndef Sudden_Underflow
  841.         if (L > 0) {
  842. #endif
  843. #ifdef IBM
  844.                 L |= Exp_msk1 >> 4;
  845. #endif
  846.                 word0(a) = L;
  847.                 word1(a) = 0;
  848. #ifndef Sudden_Underflow
  849.                 }
  850.         else {
  851.                 L = -L >> Exp_shift;
  852.                 if (L < Exp_shift) {
  853.                         word0(a) = 0x80000 >> L;
  854.                         word1(a) = 0;
  855.                         }
  856.                 else {
  857.                         word0(a) = 0;
  858.                         L -= Exp_shift;
  859.                         word1(a) = L >= 31 ? 1 : 1 << 31 - L;
  860.                         }
  861.                 }
  862. #endif
  863.         return a;
  864.         }
  865.  
  866.  static double
  867. b2d
  868. #ifdef KR_headers
  869.         (a, e) Bigint *a; int *e;
  870. #else
  871.         (Bigint *a, int *e)
  872. #endif
  873. {
  874.         unsigned long *xa, *xa0, w, y, z;
  875.         int k;
  876.         double d;
  877. #ifdef VAX
  878.         unsigned long d0, d1;
  879. #else
  880. #define d0 word0(d)
  881. #define d1 word1(d)
  882. #endif
  883.  
  884.         xa0 = a->x;
  885.         xa = xa0 + a->wds;
  886.         y = *--xa;
  887. #ifdef DEBUG
  888.         if (!y) Bug("zero y in b2d");
  889. #endif
  890.         k = hi0bits(y);
  891.         *e = 32 - k;
  892. #ifdef Pack_32
  893.         if (k < Ebits) {
  894.                 d0 = Exp_1 | y >> Ebits - k;
  895.                 w = xa > xa0 ? *--xa : 0;
  896.                 d1 = y << (32-Ebits) + k | w >> Ebits - k;
  897.                 goto ret_d;
  898.                 }
  899.         z = xa > xa0 ? *--xa : 0;
  900.         if (k -= Ebits) {
  901.                 d0 = Exp_1 | y << k | z >> 32 - k;
  902.                 y = xa > xa0 ? *--xa : 0;
  903.                 d1 = z << k | y >> 32 - k;
  904.                 }
  905.         else {
  906.                 d0 = Exp_1 | y;
  907.                 d1 = z;
  908.                 }
  909. #else
  910.         if (k < Ebits + 16) {
  911.                 z = xa > xa0 ? *--xa : 0;
  912.                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  913.                 w = xa > xa0 ? *--xa : 0;
  914.                 y = xa > xa0 ? *--xa : 0;
  915.                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  916.                 goto ret_d;
  917.                 }
  918.         z = xa > xa0 ? *--xa : 0;
  919.         w = xa > xa0 ? *--xa : 0;
  920.         k -= Ebits + 16;
  921.         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  922.         y = xa > xa0 ? *--xa : 0;
  923.         d1 = w << k + 16 | y << k;
  924. #endif
  925.  ret_d:
  926. #ifdef VAX
  927.         word0(d) = d0 >> 16 | d0 << 16;
  928.         word1(d) = d1 >> 16 | d1 << 16;
  929. #else
  930. #undef d0
  931. #undef d1
  932. #endif
  933.         return d;
  934.         }
  935.  
  936.  static Bigint *
  937. d2b
  938. #ifdef KR_headers
  939.         (d, e, bits) double d; int *e, *bits;
  940. #else
  941.         (double d, int *e, int *bits)
  942. #endif
  943. {
  944.         Bigint *b;
  945.         int de, i, k;
  946.         unsigned long *x, y, z;
  947. #ifdef VAX
  948.         unsigned long d0, d1;
  949.         d0 = word0(d) >> 16 | word0(d) << 16;
  950.         d1 = word1(d) >> 16 | word1(d) << 16;
  951. #else
  952. #define d0 word0(d)
  953. #define d1 word1(d)
  954. #endif
  955.  
  956. #ifdef Pack_32
  957.         b = Balloc(1);
  958. #else
  959.         b = Balloc(2);
  960. #endif
  961.         x = b->x;
  962.  
  963.         z = d0 & Frac_mask;
  964.         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
  965. #ifdef Sudden_Underflow
  966.         de = (int)(d0 >> Exp_shift);
  967. #ifndef IBM
  968.         z |= Exp_msk11;
  969. #endif
  970. #else
  971.         if (de = (int)(d0 >> Exp_shift))
  972.                 z |= Exp_msk1;
  973. #endif
  974. #ifdef Pack_32
  975.         if (y = d1) {
  976.                 if (k = lo0bits(&y)) {
  977.                         x[0] = y | z << 32 - k;
  978.                         z >>= k;
  979.                         }
  980.                 else
  981.                         x[0] = y;
  982.                 i = b->wds = (x[1] = z) ? 2 : 1;
  983.                 }
  984.         else {
  985. #ifdef DEBUG
  986.                 if (!z)
  987.                         Bug("Zero passed to d2b");
  988. #endif
  989.                 k = lo0bits(&z);
  990.                 x[0] = z;
  991.                 i = b->wds = 1;
  992.                 k += 32;
  993.                 }
  994. #else
  995.         if (y = d1) {
  996.                 if (k = lo0bits(&y))
  997.                         if (k >= 16) {
  998.                                 x[0] = y | z << 32 - k & 0xffff;
  999.                                 x[1] = z >> k - 16 & 0xffff;
  1000.                                 x[2] = z >> k;
  1001.                                 i = 2;
  1002.                                 }
  1003.                         else {
  1004.                                 x[0] = y & 0xffff;
  1005.                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
  1006.                                 x[2] = z >> k & 0xffff;
  1007.                                 x[3] = z >> k+16;
  1008.                                 i = 3;
  1009.                                 }
  1010.                 else {
  1011.                         x[0] = y & 0xffff;
  1012.                         x[1] = y >> 16;
  1013.                         x[2] = z & 0xffff;
  1014.                         x[3] = z >> 16;
  1015.                         i = 3;
  1016.                         }
  1017.                 }
  1018.         else {
  1019. #ifdef DEBUG
  1020.                 if (!z)
  1021.                         Bug("Zero passed to d2b");
  1022. #endif
  1023.                 k = lo0bits(&z);
  1024.                 if (k >= 16) {
  1025.                         x[0] = z;
  1026.                         i = 0;
  1027.                         }
  1028.                 else {
  1029.                         x[0] = z & 0xffff;
  1030.                         x[1] = z >> 16;
  1031.                         i = 1;
  1032.                         }
  1033.                 k += 32;
  1034.                 }
  1035.         while(!x[i])
  1036.                 --i;
  1037.         b->wds = i + 1;
  1038. #endif
  1039. #ifndef Sudden_Underflow
  1040.         if (de) {
  1041. #endif
  1042. #ifdef IBM
  1043.                 *e = (de - Bias - (P-1) << 2) + k;
  1044.                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
  1045. #else
  1046.                 *e = de - Bias - (P-1) + k;
  1047.                 *bits = P - k;
  1048. #endif
  1049. #ifndef Sudden_Underflow
  1050.                 }
  1051.         else {
  1052.                 *e = de - Bias - (P-1) + 1 + k;
  1053. #ifdef Pack_32
  1054.                 *bits = 32*i - hi0bits(x[i-1]);
  1055. #else
  1056.                 *bits = (i+2)*16 - hi0bits(x[i]);
  1057. #endif
  1058.                 }
  1059. #endif
  1060.         return b;
  1061.         }
  1062. #undef d0
  1063. #undef d1
  1064.  
  1065.  static double
  1066. ratio
  1067. #ifdef KR_headers
  1068.         (a, b) Bigint *a, *b;
  1069. #else
  1070.         (Bigint *a, Bigint *b)
  1071. #endif
  1072. {
  1073.         double da, db;
  1074.         int k, ka, kb;
  1075.  
  1076.         da = b2d(a, &ka);
  1077.         db = b2d(b, &kb);
  1078. #ifdef Pack_32
  1079.         k = ka - kb + 32*(a->wds - b->wds);
  1080. #else
  1081.         k = ka - kb + 16*(a->wds - b->wds);
  1082. #endif
  1083. #ifdef IBM
  1084.         if (k > 0) {
  1085.                 word0(da) += (k >> 2)*Exp_msk1;
  1086.                 if (k &= 3)
  1087.                         da *= 1 << k;
  1088.                 }
  1089.         else {
  1090.                 k = -k;
  1091.                 word0(db) += (k >> 2)*Exp_msk1;
  1092.                 if (k &= 3)
  1093.                         db *= 1 << k;
  1094.                 }
  1095. #else
  1096.         if (k > 0)
  1097.                 word0(da) += k*Exp_msk1;
  1098.         else {
  1099.                 k = -k;
  1100.                 word0(db) += k*Exp_msk1;
  1101.                 }
  1102. #endif
  1103.         return da / db;
  1104.         }
  1105.  
  1106.  static double
  1107. tens[] = {
  1108.                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  1109.                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  1110.                 1e20, 1e21, 1e22
  1111. #ifdef VAX
  1112.                 , 1e23, 1e24
  1113. #endif
  1114.                 };
  1115.  
  1116.  static double
  1117. #ifdef IEEE_Arith
  1118. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  1119. static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
  1120. #define n_bigtens 5
  1121. #else
  1122. #ifdef IBM
  1123. bigtens[] = { 1e16, 1e32, 1e64 };
  1124. static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  1125. #define n_bigtens 3
  1126. #else
  1127. bigtens[] = { 1e16, 1e32 };
  1128. static double tinytens[] = { 1e-16, 1e-32 };
  1129. #define n_bigtens 2
  1130. #endif
  1131. #endif
  1132.  
  1133.  double
  1134. _Xstrtod
  1135. #ifdef KR_headers
  1136.         (s00, se) CONST char *s00; char **se;
  1137. #else
  1138.         (CONST char *s00, char **se)
  1139. #endif
  1140. {
  1141.         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
  1142.                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  1143.         CONST char *s, *s0, *s1;
  1144.         double aadj, aadj1, adj, rv, rv0;
  1145.         long L;
  1146.         unsigned long y, z;
  1147.         Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  1148.     TEST_ENDIANNESS;
  1149.         sign = nz0 = nz = 0;
  1150.         rv = 0.;
  1151.         for(s = s00;;s++) switch(*s) {
  1152.                 case '-':
  1153.                         sign = 1;
  1154.                         /* no break */
  1155.                 case '+':
  1156.                         if (*++s)
  1157.                                 goto break2;
  1158.                         /* no break */
  1159.                 case 0:
  1160.                         goto ret;
  1161.                 case '\t':
  1162.                 case '\n':
  1163.                 case '\v':
  1164.                 case '\f':
  1165.                 case '\r':
  1166.                 case ' ':
  1167.                         continue;
  1168.                 default:
  1169.                         goto break2;
  1170.                 }
  1171.  break2:
  1172.         if (*s == '0') {
  1173.                 nz0 = 1;
  1174.                 while(*++s == '0') ;
  1175.                 if (!*s)
  1176.                         goto ret;
  1177.                 }
  1178.         s0 = s;
  1179.         y = z = 0;
  1180.         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  1181.                 if (nd < 9)
  1182.                         y = 10*y + c - '0';
  1183.                 else if (nd < 16)
  1184.                         z = 10*z + c - '0';
  1185.         nd0 = nd;
  1186.         if (c == '.') {
  1187.                 c = *++s;
  1188.                 if (!nd) {
  1189.                         for(; c == '0'; c = *++s)
  1190.                                 nz++;
  1191.                         if (c > '0' && c <= '9') {
  1192.                                 s0 = s;
  1193.                                 nf += nz;
  1194.                                 nz = 0;
  1195.                                 goto have_dig;
  1196.                                 }
  1197.                         goto dig_done;
  1198.                         }
  1199.                 for(; c >= '0' && c <= '9'; c = *++s) {
  1200.  have_dig:
  1201.                         nz++;
  1202.                         if (c -= '0') {
  1203.                                 nf += nz;
  1204.                                 for(i = 1; i < nz; i++)
  1205.                                         if (nd++ < 9)
  1206.                                                 y *= 10;
  1207.                                         else if (nd <= DBL_DIG + 1)
  1208.                                                 z *= 10;
  1209.                                 if (nd++ < 9)
  1210.                                         y = 10*y + c;
  1211.                                 else if (nd <= DBL_DIG + 1)
  1212.                                         z = 10*z + c;
  1213.                                 nz = 0;
  1214.                                 }
  1215.                         }
  1216.                 }
  1217.  dig_done:
  1218.         e = 0;
  1219.         if (c == 'e' || c == 'E') {
  1220.                 if (!nd && !nz && !nz0) {
  1221.                         s = s00;
  1222.                         goto ret;
  1223.                         }
  1224.                 s00 = s;
  1225.                 esign = 0;
  1226.                 switch(c = *++s) {
  1227.                         case '-':
  1228.                                 esign = 1;
  1229.                         case '+':
  1230.                                 c = *++s;
  1231.                         }
  1232.                 if (c >= '0' && c <= '9') {
  1233.                         while(c == '0')
  1234.                                 c = *++s;
  1235.                         if (c > '0' && c <= '9') {
  1236.                                 e = c - '0';
  1237.                                 s1 = s;
  1238.                                 while((c = *++s) >= '0' && c <= '9')
  1239.                                         e = 10*e + c - '0';
  1240.                                 if (s - s1 > 8)
  1241.                                         /* Avoid confusion from exponents
  1242.                                          * so large that e might overflow.
  1243.                                          */
  1244.                                         e = 9999999;
  1245.                                 if (esign)
  1246.                                         e = -e;
  1247.                                 }
  1248.                         else
  1249.                                 e = 0;
  1250.                         }
  1251.                 else
  1252.                         s = s00;
  1253.                 }
  1254.         if (!nd) {
  1255.                 if (!nz && !nz0)
  1256.                         s = s00;
  1257.                 goto ret;
  1258.                 }
  1259.         e1 = e -= nf;
  1260.  
  1261.         /* Now we have nd0 digits, starting at s0, followed by a
  1262.          * decimal point, followed by nd-nd0 digits.  The number we're
  1263.          * after is the integer represented by those digits times
  1264.          * 10**e */
  1265.  
  1266.         if (!nd0)
  1267.                 nd0 = nd;
  1268.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  1269.         rv = y;
  1270.         if (k > 9)
  1271.                 rv = tens[k - 9] * rv + z;
  1272.         if (nd <= DBL_DIG
  1273. #ifndef RND_PRODQUOT
  1274.                 && FLT_ROUNDS == 1
  1275. #endif
  1276.                         ) {
  1277.                 if (!e)
  1278.                         goto ret;
  1279.                 if (e > 0) {
  1280.                         if (e <= Ten_pmax) {
  1281. #ifdef VAX
  1282.                                 goto vax_ovfl_check;
  1283. #else
  1284.                                 /* rv = */ rounded_product(rv, tens[e]);
  1285.                                 goto ret;
  1286. #endif
  1287.                                 }
  1288.                         i = DBL_DIG - nd;
  1289.                         if (e <= Ten_pmax + i) {
  1290.                                 /* A fancier test would sometimes let us do
  1291.                                  * this for larger i values.
  1292.                                  */
  1293.                                 e -= i;
  1294.                                 rv *= tens[i];
  1295. #ifdef VAX
  1296.                                 /* VAX exponent range is so narrow we must
  1297.                                  * worry about overflow here...
  1298.                                  */
  1299.  vax_ovfl_check:
  1300.                                 word0(rv) -= P*Exp_msk1;
  1301.                                 /* rv = */ rounded_product(rv, tens[e]);
  1302.                                 if ((word0(rv) & Exp_mask)
  1303.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  1304.                                         goto ovfl;
  1305.                                 word0(rv) += P*Exp_msk1;
  1306. #else
  1307.                                 /* rv = */ rounded_product(rv, tens[e]);
  1308. #endif
  1309.                                 goto ret;
  1310.                                 }
  1311.                         }
  1312. #ifndef Inaccurate_Divide
  1313.                 else if (e >= -Ten_pmax) {
  1314.                         /* rv = */ rounded_quotient(rv, tens[-e]);
  1315.                         goto ret;
  1316.                         }
  1317. #endif
  1318.                 }
  1319.         e1 += nd - k;
  1320.  
  1321.         /* Get starting approximation = rv * 10**e1 */
  1322.  
  1323.         if (e1 > 0) {
  1324.                 if (i = e1 & 15)
  1325.                         rv *= tens[i];
  1326.                 if (e1 &= ~15) {
  1327.                         if (e1 > DBL_MAX_10_EXP) {
  1328.  ovfl:
  1329.                                 errno = ERANGE;
  1330. #ifndef HUGE_VAL
  1331. #define HUGE_VAL        1.7976931348623157E+308
  1332. #endif
  1333.                                 rv = HUGE_VAL;
  1334.                                 goto ret;
  1335.                                 }
  1336.                         if (e1 >>= 4) {
  1337.                                 for(j = 0; e1 > 1; j++, e1 >>= 1)
  1338.                                         if (e1 & 1)
  1339.                                                 rv *= bigtens[j];
  1340.                         /* The last multiplication could overflow. */
  1341.                                 word0(rv) -= P*Exp_msk1;
  1342.                                 rv *= bigtens[j];
  1343.                                 if ((z = word0(rv) & Exp_mask)
  1344.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  1345.                                         goto ovfl;
  1346.                                 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  1347.                                         /* set to largest number */
  1348.                                         /* (Can't trust DBL_MAX) */
  1349.                                         word0(rv) = Big0;
  1350.                                         word1(rv) = Big1;
  1351.                                         }
  1352.                                 else
  1353.                                         word0(rv) += P*Exp_msk1;
  1354.                                 }
  1355.  
  1356.                         }
  1357.                 }
  1358.         else if (e1 < 0) {
  1359.                 e1 = -e1;
  1360.                 if (i = e1 & 15)
  1361.                         rv /= tens[i];
  1362.                 if (e1 &= ~15) {
  1363.                         e1 >>= 4;
  1364.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  1365.                                 if (e1 & 1)
  1366.                                         rv *= tinytens[j];
  1367.                         /* The last multiplication could underflow. */
  1368.                         rv0 = rv;
  1369.                         rv *= tinytens[j];
  1370.                         if (!rv) {
  1371.                                 rv = 2.*rv0;
  1372.                                 rv *= tinytens[j];
  1373.                                 if (!rv) {
  1374.  undfl:
  1375.                                         rv = 0.;
  1376.                                         errno = ERANGE;
  1377.                                         goto ret;
  1378.                                         }
  1379.                                 word0(rv) = Tiny0;
  1380.                                 word1(rv) = Tiny1;
  1381.                                 /* The refinement below will clean
  1382.                                  * this approximation up.
  1383.                                  */
  1384.                                 }
  1385.                         }
  1386.                 }
  1387.  
  1388.         /* Now the hard part -- adjusting rv to the correct value.*/
  1389.  
  1390.         /* Put digits into bd: true value = bd * 10^e */
  1391.  
  1392.         bd0 = s2b(s0, nd0, nd, y);
  1393.  
  1394.         for(;;) {
  1395.                 bd = Balloc(bd0->k);
  1396.                 Bcopy(bd, bd0);
  1397.                 bb = d2b(rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
  1398.                 bs = i2b(1);
  1399.  
  1400.                 if (e >= 0) {
  1401.                         bb2 = bb5 = 0;
  1402.                         bd2 = bd5 = e;
  1403.                         }
  1404.                 else {
  1405.                         bb2 = bb5 = -e;
  1406.                         bd2 = bd5 = 0;
  1407.                         }
  1408.                 if (bbe >= 0)
  1409.                         bb2 += bbe;
  1410.                 else
  1411.                         bd2 -= bbe;
  1412.                 bs2 = bb2;
  1413. #ifdef Sudden_Underflow
  1414. #ifdef IBM
  1415.                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  1416. #else
  1417.                 j = P + 1 - bbbits;
  1418. #endif
  1419. #else
  1420.                 i = bbe + bbbits - 1;   /* logb(rv) */
  1421.                 if (i < Emin)   /* denormal */
  1422.                         j = bbe + (P-Emin);
  1423.                 else
  1424.                         j = P + 1 - bbbits;
  1425. #endif
  1426.                 bb2 += j;
  1427.                 bd2 += j;
  1428.                 i = bb2 < bd2 ? bb2 : bd2;
  1429.                 if (i > bs2)
  1430.                         i = bs2;
  1431.                 if (i > 0) {
  1432.                         bb2 -= i;
  1433.                         bd2 -= i;
  1434.                         bs2 -= i;
  1435.                         }
  1436.                 if (bb5 > 0) {
  1437.                         bs = pow5mult(bs, bb5);
  1438.                         bb1 = mult(bs, bb);
  1439.                         Bfree(bb);
  1440.                         bb = bb1;
  1441.                         }
  1442.                 if (bb2 > 0)
  1443.                         bb = lshift(bb, bb2);
  1444.                 if (bd5 > 0)
  1445.                         bd = pow5mult(bd, bd5);
  1446.                 if (bd2 > 0)
  1447.                         bd = lshift(bd, bd2);
  1448.                 if (bs2 > 0)
  1449.                         bs = lshift(bs, bs2);
  1450.                 delta = diff(bb, bd);
  1451.                 dsign = delta->sign;
  1452.                 delta->sign = 0;
  1453.                 i = cmp(delta, bs);
  1454.                 if (i < 0) {
  1455.                         /* Error is less than half an ulp -- check for
  1456.                          * special case of mantissa a power of two.
  1457.                          */
  1458.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask)
  1459.                                 break;
  1460.                         delta = lshift(delta,Log2P);
  1461.                         if (cmp(delta, bs) > 0)
  1462.                                 goto drop_down;
  1463.                         break;
  1464.                         }
  1465.                 if (i == 0) {
  1466.                         /* exactly half-way between */
  1467.                         if (dsign) {
  1468.                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
  1469.                                  &&  word1(rv) == 0xffffffff) {
  1470.                                         /*boundary case -- increment exponent*/
  1471.                                         word0(rv) = (word0(rv) & Exp_mask)
  1472.                                                 + Exp_msk1
  1473. #ifdef IBM
  1474.                                                 | Exp_msk1 >> 4
  1475. #endif
  1476.                                                 ;
  1477.                                         word1(rv) = 0;
  1478.                                         break;
  1479.                                         }
  1480.                                 }
  1481.                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
  1482.  drop_down:
  1483.                                 /* boundary case -- decrement exponent */
  1484. #ifdef Sudden_Underflow
  1485.                                 L = word0(rv) & Exp_mask;
  1486. #ifdef IBM
  1487.                                 if (L <  Exp_msk1)
  1488. #else
  1489.                                 if (L <= Exp_msk1)
  1490. #endif
  1491.                                         goto undfl;
  1492.                                 L -= Exp_msk1;
  1493. #else
  1494.                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
  1495. #endif
  1496.                                 word0(rv) = L | Bndry_mask1;
  1497.                                 word1(rv) = 0xffffffff;
  1498. #ifdef IBM
  1499.                                 goto cont;
  1500. #else
  1501.                                 break;
  1502. #endif
  1503.                                 }
  1504. #ifndef ROUND_BIASED
  1505.                         if (!(word1(rv) & LSB))
  1506.                                 break;
  1507. #endif
  1508.                         if (dsign)
  1509.                                 rv += ulp(rv);
  1510. #ifndef ROUND_BIASED
  1511.                         else {
  1512.                                 rv -= ulp(rv);
  1513. #ifndef Sudden_Underflow
  1514.                                 if (!rv)
  1515.                                         goto undfl;
  1516. #endif
  1517.                                 }
  1518. #endif
  1519.                         break;
  1520.                         }
  1521.                 if ((aadj = ratio(delta, bs)) <= 2.) {
  1522.                         if (dsign)
  1523.                                 aadj = aadj1 = 1.;
  1524.                         else if (word1(rv) || word0(rv) & Bndry_mask) {
  1525. #ifndef Sudden_Underflow
  1526.                                 if (word1(rv) == Tiny1 && !word0(rv))
  1527.                                         goto undfl;
  1528. #endif
  1529.                                 aadj = 1.;
  1530.                                 aadj1 = -1.;
  1531.                                 }
  1532.                         else {
  1533.                                 /* special case -- power of FLT_RADIX to be */
  1534.                                 /* rounded down... */
  1535.  
  1536.                                 if (aadj < 2./FLT_RADIX)
  1537.                                         aadj = 1./FLT_RADIX;
  1538.                                 else
  1539.                                         aadj *= 0.5;
  1540.                                 aadj1 = -aadj;
  1541.                                 }
  1542.                         }
  1543.                 else {
  1544.                         aadj *= 0.5;
  1545.                         aadj1 = dsign ? aadj : -aadj;
  1546. #ifdef Check_FLT_ROUNDS
  1547.                         switch(FLT_ROUNDS) {
  1548.                                 case 2: /* towards +infinity */
  1549.                                         aadj1 -= 0.5;
  1550.                                         break;
  1551.                                 case 0: /* towards 0 */
  1552.                                 case 3: /* towards -infinity */
  1553.                                         aadj1 += 0.5;
  1554.                                 }
  1555. #else
  1556.                         if (FLT_ROUNDS == 0)
  1557.                                 aadj1 += 0.5;
  1558. #endif
  1559.                         }
  1560.                 y = word0(rv) & Exp_mask;
  1561.  
  1562.                 /* Check for overflow */
  1563.  
  1564.                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1565.                         rv0 = rv;
  1566.                         word0(rv) -= P*Exp_msk1;
  1567.                         adj = aadj1 * ulp(rv);
  1568.                         rv += adj;
  1569.                         if ((word0(rv) & Exp_mask) >=
  1570.                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1571.                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
  1572.                                         goto ovfl;
  1573.                                 word0(rv) = Big0;
  1574.                                 word1(rv) = Big1;
  1575.                                 goto cont;
  1576.                                 }
  1577.                         else
  1578.                                 word0(rv) += P*Exp_msk1;
  1579.                         }
  1580.                 else {
  1581. #ifdef Sudden_Underflow
  1582.                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
  1583.                                 rv0 = rv;
  1584.                                 word0(rv) += P*Exp_msk1;
  1585.                                 adj = aadj1 * ulp(rv);
  1586.                                 rv += adj;
  1587. #ifdef IBM
  1588.                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
  1589. #else
  1590.                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
  1591. #endif
  1592.                                         {
  1593.                                         if (word0(rv0) == Tiny0
  1594.                                          && word1(rv0) == Tiny1)
  1595.                                                 goto undfl;
  1596.                                         word0(rv) = Tiny0;
  1597.                                         word1(rv) = Tiny1;
  1598.                                         goto cont;
  1599.                                         }
  1600.                                 else
  1601.                                         word0(rv) -= P*Exp_msk1;
  1602.                                 }
  1603.                         else {
  1604.                                 adj = aadj1 * ulp(rv);
  1605.                                 rv += adj;
  1606.                                 }
  1607. #else
  1608.                         /* Compute adj so that the IEEE rounding rules will
  1609.                          * correctly round rv + adj in some half-way cases.
  1610.                          * If rv * ulp(rv) is denormalized (i.e.,
  1611.                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1612.                          * trouble from bits lost to denormalization;
  1613.                          * example: 1.2e-307 .
  1614.                          */
  1615.                         if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
  1616.                                 aadj1 = (double)(int)(aadj + 0.5);
  1617.                                 if (!dsign)
  1618.                                         aadj1 = -aadj1;
  1619.                                 }
  1620.                         adj = aadj1 * ulp(rv);
  1621.                         rv += adj;
  1622. #endif
  1623.                         }
  1624.                 z = word0(rv) & Exp_mask;
  1625.                 if (y == z) {
  1626.                         /* Can we stop now? */
  1627.                         L = (long)aadj;
  1628.                         aadj -= L;
  1629.                         /* The tolerances below are conservative. */
  1630.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
  1631.                                 if (aadj < .4999999 || aadj > .5000001)
  1632.                                         break;
  1633.                                 }
  1634.                         else if (aadj < .4999999/FLT_RADIX)
  1635.                                 break;
  1636.                         }
  1637.  cont:
  1638.                 Bfree(bb);
  1639.                 Bfree(bd);
  1640.                 Bfree(bs);
  1641.                 Bfree(delta);
  1642.                 }
  1643.         Bfree(bb);
  1644.         Bfree(bd);
  1645.         Bfree(bs);
  1646.         Bfree(bd0);
  1647.         Bfree(delta);
  1648.  ret:
  1649.         if (se)
  1650.                 *se = (char *)s;
  1651.         return sign ? -rv : rv;
  1652.         }
  1653.  
  1654.  static int
  1655. quorem
  1656. #ifdef KR_headers
  1657.         (b, S) Bigint *b, *S;
  1658. #else
  1659.         (Bigint *b, Bigint *S)
  1660. #endif
  1661. {
  1662.         int n;
  1663.         long borrow, y;
  1664.         unsigned long carry, q, ys;
  1665.         unsigned long *bx, *bxe, *sx, *sxe;
  1666. #ifdef Pack_32
  1667.         long z;
  1668.         unsigned long si, zs;
  1669. #endif
  1670.  
  1671.         n = S->wds;
  1672. #ifdef DEBUG
  1673.         /*debug*/ if (b->wds > n)
  1674.         /*debug*/       Bug("oversize b in quorem");
  1675. #endif
  1676.         if (b->wds < n)
  1677.                 return 0;
  1678.         sx = S->x;
  1679.         sxe = sx + --n;
  1680.         bx = b->x;
  1681.         bxe = bx + n;
  1682.         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
  1683. #ifdef DEBUG
  1684.         /*debug*/ if (q > 9)
  1685.         /*debug*/       Bug("oversized quotient in quorem");
  1686. #endif
  1687.         if (q) {
  1688.                 borrow = 0;
  1689.                 carry = 0;
  1690.                 do {
  1691. #ifdef Pack_32
  1692.                         si = *sx++;
  1693.                         ys = (si & 0xffff) * q + carry;
  1694.                         zs = (si >> 16) * q + (ys >> 16);
  1695.                         carry = zs >> 16;
  1696.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1697.                         borrow = y >> 16;
  1698.                         Sign_Extend(borrow, y);
  1699.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1700.                         borrow = z >> 16;
  1701.                         Sign_Extend(borrow, z);
  1702.                         Storeinc(bx, z, y);
  1703. #else
  1704.                         ys = *sx++ * q + carry;
  1705.                         carry = ys >> 16;
  1706.                         y = *bx - (ys & 0xffff) + borrow;
  1707.                         borrow = y >> 16;
  1708.                         Sign_Extend(borrow, y);
  1709.                         *bx++ = y & 0xffff;
  1710. #endif
  1711.                         }
  1712.                         while(sx <= sxe);
  1713.                 if (!*bxe) {
  1714.                         bx = b->x;
  1715.                         while(--bxe > bx && !*bxe)
  1716.                                 --n;
  1717.                         b->wds = n;
  1718.                         }
  1719.                 }
  1720.         if (cmp(b, S) >= 0) {
  1721.                 q++;
  1722.                 borrow = 0;
  1723.                 carry = 0;
  1724.                 bx = b->x;
  1725.                 sx = S->x;
  1726.                 do {
  1727. #ifdef Pack_32
  1728.                         si = *sx++;
  1729.                         ys = (si & 0xffff) + carry;
  1730.                         zs = (si >> 16) + (ys >> 16);
  1731.                         carry = zs >> 16;
  1732.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1733.                         borrow = y >> 16;
  1734.                         Sign_Extend(borrow, y);
  1735.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1736.                         borrow = z >> 16;
  1737.                         Sign_Extend(borrow, z);
  1738.                         Storeinc(bx, z, y);
  1739. #else
  1740.                         ys = *sx++ + carry;
  1741.                         carry = ys >> 16;
  1742.                         y = *bx - (ys & 0xffff) + borrow;
  1743.                         borrow = y >> 16;
  1744.                         Sign_Extend(borrow, y);
  1745.                         *bx++ = y & 0xffff;
  1746. #endif
  1747.                         }
  1748.                         while(sx <= sxe);
  1749.                 bx = b->x;
  1750.                 bxe = bx + n;
  1751.                 if (!*bxe) {
  1752.                         while(--bxe > bx && !*bxe)
  1753.                                 --n;
  1754.                         b->wds = n;
  1755.                         }
  1756.                 }
  1757.         return q;
  1758.         }
  1759.  
  1760. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  1761.  *
  1762.  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  1763.  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
  1764.  *
  1765.  * Modifications:
  1766.  *      1. Rather than iterating, we use a simple numeric overestimate
  1767.  *         to determine k = floor(log10(d)).  We scale relevant
  1768.  *         quantities using O(log2(k)) rather than O(k) multiplications.
  1769.  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  1770.  *         try to generate digits strictly left to right.  Instead, we
  1771.  *         compute with fewer bits and propagate the carry if necessary
  1772.  *         when rounding the final digit up.  This is often faster.
  1773.  *      3. Under the assumption that input will be rounded nearest,
  1774.  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  1775.  *         That is, we allow equality in stopping tests when the
  1776.  *         round-nearest rule will give the same floating-point value
  1777.  *         as would satisfaction of the stopping test with strict
  1778.  *         inequality.
  1779.  *      4. We remove common factors of powers of 2 from relevant
  1780.  *         quantities.
  1781.  *      5. When converting floating-point integers less than 1e16,
  1782.  *         we use floating-point arithmetic rather than resorting
  1783.  *         to multiple-precision integers.
  1784.  *      6. When asked to produce fewer than 15 digits, we first try
  1785.  *         to get by with floating-point arithmetic; we resort to
  1786.  *         multiple-precision integer arithmetic only if we cannot
  1787.  *         guarantee that the floating-point calculation has given
  1788.  *         the correctly rounded result.  For k requested digits and
  1789.  *         "uniformly" distributed input, the probability is
  1790.  *         something like 10^(k-15) that we must resort to the long
  1791.  *         calculation.
  1792.  */
  1793.  
  1794.  char *
  1795. dtoa
  1796. #ifdef KR_headers
  1797.         (d, mode, ndigits, decpt, sign, rve)
  1798.         double d; int mode, ndigits, *decpt, *sign; char **rve;
  1799. #else
  1800.         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  1801. #endif
  1802. {
  1803.  /*     Arguments ndigits, decpt, sign are similar to those
  1804.         of ecvt and fcvt; trailing zeros are suppressed from
  1805.         the returned string.  If not null, *rve is set to point
  1806.         to the end of the return value.  If d is +-Infinity or NaN,
  1807.         then *decpt is set to 9999.
  1808.  
  1809.         mode:
  1810.                 0 ==> shortest string that yields d when read in
  1811.                         and rounded to nearest.
  1812.                 1 ==> like 0, but with Steele & White stopping rule;
  1813.                         e.g. with IEEE P754 arithmetic , mode 0 gives
  1814.                         1e23 whereas mode 1 gives 9.999999999999999e22.
  1815.                 2 ==> max(1,ndigits) significant digits.  This gives a
  1816.                         return value similar to that of ecvt, except
  1817.                         that trailing zeros are suppressed.
  1818.                 3 ==> through ndigits past the decimal point.  This
  1819.                         gives a return value similar to that from fcvt,
  1820.                         except that trailing zeros are suppressed, and
  1821.                         ndigits can be negative.
  1822.                 4-9 should give the same return values as 2-3, i.e.,
  1823.                         4 <= mode <= 9 ==> same return as mode
  1824.                         2 + (mode & 1).  These modes are mainly for
  1825.                         debugging; often they run slower but sometimes
  1826.                         faster than modes 2-3.
  1827.                 4,5,8,9 ==> left-to-right digit generation.
  1828.                 6-9 ==> don't try fast floating-point estimate
  1829.                         (if applicable).
  1830.  
  1831.                 Values of mode other than 0-9 are treated as mode 0.
  1832.  
  1833.                 Sufficient space is allocated to the return value
  1834.                 to hold the suppressed trailing zeros.
  1835.         */
  1836.  
  1837.         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
  1838.                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
  1839.                 spec_case, try_quick;
  1840.         long L;
  1841. #ifndef Sudden_Underflow
  1842.         int denorm;
  1843.         unsigned long x;
  1844. #endif
  1845.         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
  1846.         double d2, ds, eps;
  1847.         char *s, *s0;
  1848.         static Bigint *result;
  1849.         static int result_k;
  1850.  
  1851.     TEST_ENDIANNESS;
  1852.         if (result) {
  1853.                 result->k = result_k;
  1854.                 result->maxwds = 1 << result_k;
  1855.                 Bfree(result);
  1856.                 result = 0;
  1857.                 }
  1858.  
  1859.         if (word0(d) & Sign_bit) {
  1860.                 /* set sign for everything, including 0's and NaNs */
  1861.                 *sign = 1;
  1862.                 word0(d) &= ~Sign_bit;  /* clear sign bit */
  1863.                 }
  1864.         else
  1865.                 *sign = 0;
  1866.  
  1867. #if defined(IEEE_Arith) + defined(VAX)
  1868. #ifdef IEEE_Arith
  1869.         if ((word0(d) & Exp_mask) == Exp_mask)
  1870. #else
  1871.         if (word0(d)  == 0x8000)
  1872. #endif
  1873.                 {
  1874.                 /* Infinity or NaN */
  1875.                 *decpt = 9999;
  1876.                 s =
  1877. #ifdef IEEE_Arith
  1878.                         !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
  1879. #endif
  1880.                                 "NaN";
  1881.                 if (rve)
  1882.                         *rve =
  1883. #ifdef IEEE_Arith
  1884.                                 s[3] ? s + 8 :
  1885. #endif
  1886.                                                 s + 3;
  1887.                 return s;
  1888.                 }
  1889. #endif
  1890. #ifdef IBM
  1891.         d += 0; /* normalize */
  1892. #endif
  1893.         if (!d) {
  1894.                 *decpt = 1;
  1895.                 s = "0";
  1896.                 if (rve)
  1897.                         *rve = s + 1;
  1898.                 return s;
  1899.                 }
  1900.  
  1901.         b = d2b(d, &be, &bbits);
  1902. #ifdef Sudden_Underflow
  1903.         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
  1904. #else
  1905.         if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
  1906. #endif
  1907.                 d2 = d;
  1908.                 word0(d2) &= Frac_mask1;
  1909.                 word0(d2) |= Exp_11;
  1910. #ifdef IBM
  1911.                 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
  1912.                         d2 /= 1 << j;
  1913. #endif
  1914.  
  1915.                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
  1916.                  * log10(x)      =  log(x) / log(10)
  1917.                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  1918.                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
  1919.                  *
  1920.                  * This suggests computing an approximation k to log10(d) by
  1921.                  *
  1922.                  * k = (i - Bias)*0.301029995663981
  1923.                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  1924.                  *
  1925.                  * We want k to be too large rather than too small.
  1926.                  * The error in the first-order Taylor series approximation
  1927.                  * is in our favor, so we just round up the constant enough
  1928.                  * to compensate for any error in the multiplication of
  1929.                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
  1930.                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  1931.                  * adding 1e-13 to the constant term more than suffices.
  1932.                  * Hence we adjust the constant term to 0.1760912590558.
  1933.                  * (We could get a more accurate k by invoking log10,
  1934.                  *  but this is probably not worthwhile.)
  1935.                  */
  1936.  
  1937.                 i -= Bias;
  1938. #ifdef IBM
  1939.                 i <<= 2;
  1940.                 i += j;
  1941. #endif
  1942. #ifndef Sudden_Underflow
  1943.                 denorm = 0;
  1944.                 }
  1945.         else {
  1946.                 /* d is denormalized */
  1947.  
  1948.                 i = bbits + be + (Bias + (P-1) - 1);
  1949.                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
  1950.                             : word1(d) << 32 - i;
  1951.                 d2 = x;
  1952.                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
  1953.                 i -= (Bias + (P-1) - 1) + 1;
  1954.                 denorm = 1;
  1955.                 }
  1956. #endif
  1957.         ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
  1958.         k = (int)ds;
  1959.         if (ds < 0. && ds != k)
  1960.                 k--;    /* want k = floor(ds) */
  1961.         k_check = 1;
  1962.         if (k >= 0 && k <= Ten_pmax) {
  1963.                 if (d < tens[k])
  1964.                         k--;
  1965.                 k_check = 0;
  1966.                 }
  1967.         j = bbits - i - 1;
  1968.         if (j >= 0) {
  1969.                 b2 = 0;
  1970.                 s2 = j;
  1971.                 }
  1972.         else {
  1973.                 b2 = -j;
  1974.                 s2 = 0;
  1975.                 }
  1976.         if (k >= 0) {
  1977.                 b5 = 0;
  1978.                 s5 = k;
  1979.                 s2 += k;
  1980.                 }
  1981.         else {
  1982.                 b2 -= k;
  1983.                 b5 = -k;
  1984.                 s5 = 0;
  1985.                 }
  1986.         if (mode < 0 || mode > 9)
  1987.                 mode = 0;
  1988.         try_quick = 1;
  1989.         if (mode > 5) {
  1990.                 mode -= 4;
  1991.                 try_quick = 0;
  1992.                 }
  1993.         leftright = 1;
  1994.         switch(mode) {
  1995.                 case 0:
  1996.                 case 1:
  1997.                         ilim = ilim1 = -1;
  1998.                         i = 18;
  1999.                         ndigits = 0;
  2000.                         break;
  2001.                 case 2:
  2002.                         leftright = 0;
  2003.                         /* no break */
  2004.                 case 4:
  2005.                         if (ndigits <= 0)
  2006.                                 ndigits = 1;
  2007.                         ilim = ilim1 = i = ndigits;
  2008.                         break;
  2009.                 case 3:
  2010.                         leftright = 0;
  2011.                         /* no break */
  2012.                 case 5:
  2013.                         i = ndigits + k + 1;
  2014.                         ilim = i;
  2015.                         ilim1 = i - 1;
  2016.                         if (i <= 0)
  2017.                                 i = 1;
  2018.                 }
  2019.         j = sizeof(unsigned long);
  2020.         for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j <= i;
  2021.                 j <<= 1) result_k++;
  2022.         result = Balloc(result_k);
  2023.         s = s0 = (char *)result;
  2024.  
  2025.         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
  2026.  
  2027.                 /* Try to get by with floating-point arithmetic. */
  2028.  
  2029.                 i = 0;
  2030.                 d2 = d;
  2031.                 k0 = k;
  2032.                 ilim0 = ilim;
  2033.                 ieps = 2; /* conservative */
  2034.                 if (k > 0) {
  2035.                         ds = tens[k&0xf];
  2036.                         j = k >> 4;
  2037.                         if (j & Bletch) {
  2038.                                 /* prevent overflows */
  2039.                                 j &= Bletch - 1;
  2040.                                 d /= bigtens[n_bigtens-1];
  2041.                                 ieps++;
  2042.                                 }
  2043.                         for(; j; j >>= 1, i++)
  2044.                                 if (j & 1) {
  2045.                                         ieps++;
  2046.                                         ds *= bigtens[i];
  2047.                                         }
  2048.                         d /= ds;
  2049.                         }
  2050.                 else if (j1 = -k) {
  2051.                         d *= tens[j1 & 0xf];
  2052.                         for(j = j1 >> 4; j; j >>= 1, i++)
  2053.                                 if (j & 1) {
  2054.                                         ieps++;
  2055.                                         d *= bigtens[i];
  2056.                                         }
  2057.                         }
  2058.                 if (k_check && d < 1. && ilim > 0) {
  2059.                         if (ilim1 <= 0)
  2060.                                 goto fast_failed;
  2061.                         ilim = ilim1;
  2062.                         k--;
  2063.                         d *= 10.;
  2064.                         ieps++;
  2065.                         }
  2066.                 eps = ieps*d + 7.;
  2067.                 word0(eps) -= (P-1)*Exp_msk1;
  2068.                 if (ilim == 0) {
  2069.                         S = mhi = 0;
  2070.                         d -= 5.;
  2071.                         if (d > eps)
  2072.                                 goto one_digit;
  2073.                         if (d < -eps)
  2074.                                 goto no_digits;
  2075.                         goto fast_failed;
  2076.                         }
  2077. #ifndef No_leftright
  2078.                 if (leftright) {
  2079.                         /* Use Steele & White method of only
  2080.                          * generating digits needed.
  2081.                          */
  2082.                         eps = 0.5/tens[ilim-1] - eps;
  2083.                         for(i = 0;;) {
  2084.                                 L = (long)d;
  2085.                                 d -= L;
  2086.                                 *s++ = '0' + (int)L;
  2087.                                 if (d < eps)
  2088.                                         goto ret1;
  2089.                                 if (1. - d < eps)
  2090.                                         goto bump_up;
  2091.                                 if (++i >= ilim)
  2092.                                         break;
  2093.                                 eps *= 10.;
  2094.                                 d *= 10.;
  2095.                                 }
  2096.                         }
  2097.                 else {
  2098. #endif
  2099.                         /* Generate ilim digits, then fix them up. */
  2100.                         eps *= tens[ilim-1];
  2101.                         for(i = 1;; i++, d *= 10.) {
  2102.                                 L = (long)d;
  2103.                                 d -= L;
  2104.                                 *s++ = '0' + (int)L;
  2105.                                 if (i == ilim) {
  2106.                                         if (d > 0.5 + eps)
  2107.                                                 goto bump_up;
  2108.                                         else if (d < 0.5 - eps) {
  2109.                                                 while(*--s == '0');
  2110.                                                 s++;
  2111.                                                 goto ret1;
  2112.                                                 }
  2113.                                         break;
  2114.                                         }
  2115.                                 }
  2116. #ifndef No_leftright
  2117.                         }
  2118. #endif
  2119.  fast_failed:
  2120.                 s = s0;
  2121.                 d = d2;
  2122.                 k = k0;
  2123.                 ilim = ilim0;
  2124.                 }
  2125.  
  2126.         /* Do we have a "small" integer? */
  2127.  
  2128.         if (be >= 0 && k <= Int_max) {
  2129.                 /* Yes. */
  2130.                 ds = tens[k];
  2131.                 if (ndigits < 0 && ilim <= 0) {
  2132.                         S = mhi = 0;
  2133.                         if (ilim < 0 || d <= 5*ds)
  2134.                                 goto no_digits;
  2135.                         goto one_digit;
  2136.                         }
  2137.                 for(i = 1;; i++) {
  2138.                         L = (long)(d / ds);
  2139.                         d -= L*ds;
  2140. #ifdef Check_FLT_ROUNDS
  2141.                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  2142.                         if (d < 0) {
  2143.                                 L--;
  2144.                                 d += ds;
  2145.                                 }
  2146. #endif
  2147.                         *s++ = '0' + (int)L;
  2148.                         if (i == ilim) {
  2149.                                 d += d;
  2150.                                 if (d > ds || d == ds && L & 1) {
  2151.  bump_up:
  2152.                                         while(*--s == '9')
  2153.                                                 if (s == s0) {
  2154.                                                         k++;
  2155.                                                         *s = '0';
  2156.                                                         break;
  2157.                                                         }
  2158.                                         ++*s++;
  2159.                                         }
  2160.                                 break;
  2161.                                 }
  2162.                         if (!(d *= 10.))
  2163.                                 break;
  2164.                         }
  2165.                 goto ret1;
  2166.                 }
  2167.  
  2168.         m2 = b2;
  2169.         m5 = b5;
  2170.         mhi = mlo = 0;
  2171.         if (leftright) {
  2172.                 if (mode < 2) {
  2173.                         i =
  2174. #ifndef Sudden_Underflow
  2175.                                 denorm ? be + (Bias + (P-1) - 1 + 1) :
  2176. #endif
  2177. #ifdef IBM
  2178.                                 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
  2179. #else
  2180.                                 1 + P - bbits;
  2181. #endif
  2182.                         }
  2183.                 else {
  2184.                         j = ilim - 1;
  2185.                         if (m5 >= j)
  2186.                                 m5 -= j;
  2187.                         else {
  2188.                                 s5 += j -= m5;
  2189.                                 b5 += j;
  2190.                                 m5 = 0;
  2191.                                 }
  2192.                         if ((i = ilim) < 0) {
  2193.                                 m2 -= i;
  2194.                                 i = 0;
  2195.                                 }
  2196.                         }
  2197.                 b2 += i;
  2198.                 s2 += i;
  2199.                 mhi = i2b(1);
  2200.                 }
  2201.         if (m2 > 0 && s2 > 0) {
  2202.                 i = m2 < s2 ? m2 : s2;
  2203.                 b2 -= i;
  2204.                 m2 -= i;
  2205.                 s2 -= i;
  2206.                 }
  2207.         if (b5 > 0) {
  2208.                 if (leftright) {
  2209.                         if (m5 > 0) {
  2210.                                 mhi = pow5mult(mhi, m5);
  2211.                                 b1 = mult(mhi, b);
  2212.                                 Bfree(b);
  2213.                                 b = b1;
  2214.                                 }
  2215.                         if (j = b5 - m5)
  2216.                                 b = pow5mult(b, j);
  2217.                         }
  2218.                 else
  2219.                         b = pow5mult(b, b5);
  2220.                 }
  2221.         S = i2b(1);
  2222.         if (s5 > 0)
  2223.                 S = pow5mult(S, s5);
  2224.  
  2225.         /* Check for special case that d is a normalized power of 2. */
  2226.  
  2227.         if (mode < 2) {
  2228.                 if (!word1(d) && !(word0(d) & Bndry_mask)
  2229. #ifndef Sudden_Underflow
  2230.                  && word0(d) & Exp_mask
  2231. #endif
  2232.                                 ) {
  2233.                         /* The special case */
  2234.                         b2 += Log2P;
  2235.                         s2 += Log2P;
  2236.                         spec_case = 1;
  2237.                         }
  2238.                 else
  2239.                         spec_case = 0;
  2240.                 }
  2241.  
  2242.         /* Arrange for convenient computation of quotients:
  2243.          * shift left if necessary so divisor has 4 leading 0 bits.
  2244.          *
  2245.          * Perhaps we should just compute leading 28 bits of S once
  2246.          * and for all and pass them and a shift to quorem, so it
  2247.          * can do shifts and ors to compute the numerator for q.
  2248.          */
  2249. #ifdef Pack_32
  2250.         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
  2251.                 i = 32 - i;
  2252. #else
  2253.         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
  2254.                 i = 16 - i;
  2255. #endif
  2256.         if (i > 4) {
  2257.                 i -= 4;
  2258.                 b2 += i;
  2259.                 m2 += i;
  2260.                 s2 += i;
  2261.                 }
  2262.         else if (i < 4) {
  2263.                 i += 28;
  2264.                 b2 += i;
  2265.                 m2 += i;
  2266.                 s2 += i;
  2267.                 }
  2268.         if (b2 > 0)
  2269.                 b = lshift(b, b2);
  2270.         if (s2 > 0)
  2271.                 S = lshift(S, s2);
  2272.         if (k_check) {
  2273.                 if (cmp(b,S) < 0) {
  2274.                         k--;
  2275.                         b = multadd(b, 10, 0);  /* we botched the k estimate */
  2276.                         if (leftright)
  2277.                                 mhi = multadd(mhi, 10, 0);
  2278.                         ilim = ilim1;
  2279.                         }
  2280.                 }
  2281.         if (ilim <= 0 && mode > 2) {
  2282.                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
  2283.                         /* no digits, fcvt style */
  2284.  no_digits:
  2285.                         k = -1 - ndigits;
  2286.                         goto ret;
  2287.                         }
  2288.  one_digit:
  2289.                 *s++ = '1';
  2290.                 k++;
  2291.                 goto ret;
  2292.                 }
  2293.         if (leftright) {
  2294.                 if (m2 > 0)
  2295.                         mhi = lshift(mhi, m2);
  2296.  
  2297.                 /* Compute mlo -- check for special case
  2298.                  * that d is a normalized power of 2.
  2299.                  */
  2300.  
  2301.                 mlo = mhi;
  2302.                 if (spec_case) {
  2303.                         mhi = Balloc(mhi->k);
  2304.                         Bcopy(mhi, mlo);
  2305.                         mhi = lshift(mhi, Log2P);
  2306.                         }
  2307.  
  2308.                 for(i = 1;;i++) {
  2309.                         dig = quorem(b,S) + '0';
  2310.                         /* Do we yet have the shortest decimal string
  2311.                          * that will round to d?
  2312.                          */
  2313.                         j = cmp(b, mlo);
  2314.                         delta = diff(S, mhi);
  2315.                         j1 = delta->sign ? 1 : cmp(b, delta);
  2316.                         Bfree(delta);
  2317. #ifndef ROUND_BIASED
  2318.                         if (j1 == 0 && !mode && !(word1(d) & 1)) {
  2319.                                 if (dig == '9')
  2320.                                         goto round_9_up;
  2321.                                 if (j > 0)
  2322.                                         dig++;
  2323.                                 *s++ = dig;
  2324.                                 goto ret;
  2325.                                 }
  2326. #endif
  2327.                         if (j < 0 || j == 0 && !mode
  2328. #ifndef ROUND_BIASED
  2329.                                                         && !(word1(d) & 1)
  2330. #endif
  2331.                                         ) {
  2332.                                 if (j1 > 0) {
  2333.                                         b = lshift(b, 1);
  2334.                                         j1 = cmp(b, S);
  2335.                                         if ((j1 > 0 || j1 == 0 && dig & 1)
  2336.                                         && dig++ == '9')
  2337.                                                 goto round_9_up;
  2338.                                         }
  2339.                                 *s++ = dig;
  2340.                                 goto ret;
  2341.                                 }
  2342.                         if (j1 > 0) {
  2343.                                 if (dig == '9') { /* possible if i == 1 */
  2344.  round_9_up:
  2345.                                         *s++ = '9';
  2346.                                         goto roundoff;
  2347.                                         }
  2348.                                 *s++ = dig + 1;
  2349.                                 goto ret;
  2350.                                 }
  2351.                         *s++ = dig;
  2352.                         if (i == ilim)
  2353.                                 break;
  2354.                         b = multadd(b, 10, 0);
  2355.                         if (mlo == mhi)
  2356.                                 mlo = mhi = multadd(mhi, 10, 0);
  2357.                         else {
  2358.                                 mlo = multadd(mlo, 10, 0);
  2359.                                 mhi = multadd(mhi, 10, 0);
  2360.                                 }
  2361.                         }
  2362.                 }
  2363.         else
  2364.                 for(i = 1;; i++) {
  2365.                         *s++ = dig = quorem(b,S) + '0';
  2366.                         if (i >= ilim)
  2367.                                 break;
  2368.                         b = multadd(b, 10, 0);
  2369.                         }
  2370.  
  2371.         /* Round off last digit */
  2372.  
  2373.         b = lshift(b, 1);
  2374.         j = cmp(b, S);
  2375.         if (j > 0 || j == 0 && dig & 1) {
  2376.  roundoff:
  2377.                 while(*--s == '9')
  2378.                         if (s == s0) {
  2379.                                 k++;
  2380.                                 *s++ = '1';
  2381.                                 goto ret;
  2382.                                 }
  2383.                 ++*s++;
  2384.                 }
  2385.         else {
  2386.                 while(*--s == '0');
  2387.                 s++;
  2388.                 }
  2389.  ret:
  2390.         Bfree(S);
  2391.         if (mhi) {
  2392.                 if (mlo && mlo != mhi)
  2393.                         Bfree(mlo);
  2394.                 Bfree(mhi);
  2395.                 }
  2396.  ret1:
  2397.         Bfree(b);
  2398.         *s = 0;
  2399.         *decpt = k + 1;
  2400.         if (rve)
  2401.                 *rve = s;
  2402.         return s0;
  2403.         }
  2404. #endif /* USE_DTOA */
  2405.